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ABSTRACT 

We study the phase-space behaviour of nearby trajectories in integrable potentials. 
We show that the separation of nearby orbits initially diverges very fast, mimicking a 
nearly exponential behaviour, while at late times it grows linearly. This initial expo- 
nential phase, known as Miller's instability, is commonly found in N-body simulations, 
and has been attributed to short-term (microscopic) N-body chaos. However we show 
here analytically that the initial divergence is simply due to the shape of an orbit in 
phase-space. This result confirms previous suspicions that this transient phenomenon 
is not related to an instability in the sense of non-integrable behaviour in the dynamics 
of N-body systems. 
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1 INTRODUCTION 

The problem of how exactly galaxies reach their final equi- 
librium configuration is still unsolved. It is clear that, un- 
like for gases, two-body collisions between stars in galax- 
ies are not the driving mechanism to reach a relaxed 
state, since the associated timescales are exceedingly large 
(|Binnev fc Tremaind fl987h . In an attempt to explain the 
road to equilibrium from a statistical mechanics point of 
view, Lynden-Bell (1967) introduced the concept of "vio- 
lent relaxation". In this context, the relaxation is reached 
through the effects of a "violently changing" gravita- 
tional field. However, the detai led physics of this process 
also remain to be understood l|Arad fc Lvnden-Belll 120051 ; 
IValluri et al1l2007h . 

Besides the statistical mechanics approach, it is also 
possible to study the problem of "relaxation" at the level 
of orbits. In this case, it is useful to introduce the con- 
cept of mixing, by which we mean how quickly nearby 
particle trajectories diverge in (phase-)space as a func- 
tion of time. In the case of time-independent gravita- 
tional potentials it is customary to classify mixing into 
two types. If the particles move in an integrable poten- 
tial, nearby orb it s will diverge as a power-law in time, e.g. 
iHelmi fc White] d 19991). Thi s pro cess is known as phase- 
mixing ( Binnev fc Tremaind I"l987i[). However, when the po- 
tential admits a certain amount of chaos, there exist regions 
of phase-space where nearby orbits diverge exponentially, 
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evidencing an e xtreme sensitivity to small change s in the ini- 
tial conditions (jLichtenberg fc Lieberman I 1983|). This pro- 
cess is known as chaotic-mixing i|Kandrup fc Mahonl Il994l ; 
lKandrud["l998l ). 

These mixing processes can also take place in a time- 
dependent gravitational potential, in which case the energies 
of the particles will not be constant. The degree of "sticki- 
ness" , quantified by the time-evolution of the divergence of 
nearby orbits would then measure the degree of ergodicity 
of the mixing process. In the case of chaotic mixing, this 
could lead to a system that does not have much memory of 
its evolutionary history. The timescales for evolution could 
be relatively short, and in principle, this process could be 
important in the path towards equilibrium for gal axies in 
the Universe (|Merrittll2005l : IValluri fc Merrittll2000l ). 

Since the 1970s N-body simulations have become the 
standard tool for studies of the formation and dynamics of 
structures in the Universe. The question of whether they 
are a faithful representation of the Universe has always at- 
tracted sig nificant attention. This is espec i ally t rue in re- 
cent years l|Diemand et al. 1 12004| ; iBinnev 1 12004 ) , particu- 
larly with the findi ng that dark-matter halos have univer- 
sal density profi les l|Navarro et al.lll996l : iMoore et aTlll999l ; 
IWeinberdhoOlal lbh. 

One of the first works t o focus on how N-body sys- 
tems evolve was iMilled ( 1964). Using what must have been 
the very first computers in the world, he simulated a self- 
consistent system in virial equilibrium of 8 upto 32 parti- 
cles distributed randomly in a cubic volume. Miller found 
that the trajectories of neighbouring particles initially di- 
verged exponentially. This initial transient has been con- 
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firmed using numerical e xperiments with a significantly 
larger number of particles ([L ecar 1968; Kandrup fc Smith! 



I1991I : IValluri fc Merrittll2000l : IHemsendorf fc Merrittll2002h . 
as well as with various degrees of numerical softening 
l|Kandrup fc Siderisl l200ll ). This implies that the initial ex- 
ponential divergence cannot be purely attributed to the very 
grainy nature of the gravitational potential in Miller's ex- 
periments. Furthermore, even in high-resolution N-body re- 
alizations of well-behaved integrable systems such as the 
Plummer sphere, nearby orbits experi ment a phase of ex- 
pone ntial separation at very early times (|Kandrup fc Siderisl 
120031 ). This initial exponential divergence present in N-body 
simulations is now known as "Miller's instability". Under- 
standing this puzzle is the focus of this paper. 

That N-body systems would show a certain degree 
of chaoticity is not necessarily unexpected. However, it 
seems natural to expect that the larger the number of 
particles used to represent an otherwise integrable smooth 
gravitational potential, the more faithful the represen- 
tation, and hence the lesser the degree of "numerical" 
chaos (|Quinlan fc Tremaind Il992h . There is now signifi- 
cant evidence that when such a system is represented by 
a sufficiently large number of particles, it does tend to 
the behav iour expected from the collisio n less Bolt z mann 
equation fCoodman. Heggie fc Hutl Il993l ; |E1-Zant I |2002| ; 
iKandrup fc Siderisll2003l ; ISiderisll2004l ). 

Nevertheless, even in these high-resolution exper- 
iments the initial ex ponen tial growth phase is p resent 



llKandrup fc Smith I Il99ll ; IValluri fc Merrittl 120001: 
Kan drup fc Siderisl boOll). Furthe r more, there is evidence 
ilGoodman, Heggie fc Hutl Il993l ; IHemsendorf fc Merrittl 
2002) that the rate of divergence associated to this phase 
increases in proportion to the number of particles used. 
Because Miller's instability only lasts for a very short 
timescale this does not imply that the sys tem is (macro- 
scopically) chao tic l|Valluri fc Merrittl [2000) . As stated by 
lEl-Zant 1 (|2002i ) it is likely that the "mechanism leading 
to the short e-folding time in point particle systems is 
physically unimportant" . 

So, while the existence of a continuum limit in N-body 
systems appears to be more or less established for long 
timescales, on short timescales Miller's instability remains 
a puzzle. The physical mechanism responsible for this was 
hitherto unknown. It seems quite unlikely that collisions be- 
tween particles could be important on timescales as short as 
one-tenth of th e crossing time of the system , as measured 
for example by IHemsendorf fc Mer ritt ( 2002-;) . Microscopic 
chaos arising from "white-noise" or poor orbit integrations 
are also unlikely to be important on those timescales, par- 
ticularly in integrable (well-behaved) potentials. 

In this paper, we tackle this paradox by studying the 
initial behaviour of nearby characteristics in an integrable 
smooth (and analytic) potential. Our aim is to understand 
how these nearby characteristics diverge on short timescales, 
and if they do so at nearly exponential rates. As we shall 
demonstrate below, this is indeed the case. The initial be- 
haviour mimics an exponential divergence, but since the sys- 
tem is fully integrable this is not related to the presence of 
chaos. This near-exponential behaviour merely reflects the 
time evolution of an orbit in phase-space. This result shows 
that there is no need to introduce the concept of microscopic 
chaos, and confirms previous suspicions that this transient 



phenomenon is not related to an instability in the sense of 
non-integrable behaviour in the dynamics of N-body sys- 
tems. 

In this paper we describe the evolution of nearby orbits 
in phase-space, and in particular in configuration space, ex - 
panding upon a model developed bv lHelmi fc White] (|l999l ) 
(hereafter HW). The details of this formalism are given in 
Sec. [2] In this Section we focus in detail on the behaviour of 
nearby orbits in a Plummer potential. In Sec. [3] we summa- 
rize our results. 



2 THE EVOLUTION IN PHASE-SPACE OF 
NEARBY ORBITS IN INTEGRABLE 
POTENTIALS 

The problem of the phase-space evolution of nearby orbits 
has many applications. Some of the most recent are related 
to the evolution of streams formed by the disruption of satel- 
lite systems (dwarf galaxies, globular clusters) in an external 
(Galactic) potential. This is also the basis of the formalism 
that HW developed, which is based on the conservation of 
phase-space density. It consists in a mapping from the ini- 
tial to the final configurations using adiabatic invariants (a 
schematic flow chart is given in Figure [Q. 

The basic idea is to map the initial system onto action- 
angle space, then follow the much simpler evolution in this 
space, and finally transform back locally onto observable co- 
ordinates (all these being linear transformations; for details 
see HW). This method, which uses action-angle variables, is 
very general an d can be applied to any potential that admit s 
regular orbits (iGoldsteinl Il959l ; iBinnev fc Tremainel Il987l ). 
However, if the potential is separable, the implementation is 
simpler. This includes all spherically symmetric potentials 
but only few axisymmetric and triaxi al cases, such as th e 
general class of S tackel potentials e.g. | Lvnden-Belll l|l962l ); 
|Pe Zeeuwl l| 19851 ) ; |Peionghe fc De Zeeuw I (|l988h - In this pa- 
per we shall only focus on spherical potentials because these 
are the simplest to model, while at the same time, they ev- 
idence a generic behaviour. 

Therefore, instead of following the evolution of pairs of 
nearby orbits as is traditional in N-body systems, we follow 
the evolution of a distribution function in phase-space. In 
particular, and for simplicity, we assume this distribution 
function to be a multivariate Gaussian (in 6-dimensions). 

The work presented here exploits and expands in two 
new directions the HW algorithm. Firstly, we now compute 
the behaviour of streams in physical space (to be able to 
determine the evolution of the spatial separation of nearby 
trajectories) . Secondly, we derive explicitly new analytic ex- 
pressions for this evolution on short timescales. 



2.1 The evolution of the distribution function 

As discussed above, we assume that the initial distribu- 
tion function of the system is a multivariate Gaussian in 
w — (x, v) coordinates centered on (wo) (a given particle 
or orbit): 



f{vj,to) = /oexp 



(1) 
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t = o 



Gaussian f(ra (x ,v )) 



Locally 
Gaussian f(tn(x,v),t) 



Gaussian f(w ((|) ,J o )) 



Time evolution 0(t) 

<j> = (h, + fi(J) t ■ 

J = J„ 



Gaussian f(w((|),J ,t)) 



Figure 1. Flow chart showing the basic steps of our analytic 
formalism to measure the evolution of a system in phase-space. 
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where A ro .o = w — (vjo), and cr^o is the variance matrix 
(the inverse of the covariance matrix) at the initial time: 



S x ,o 

C xv .o 



C X v,0 
O"v,0 



(2) 



For example, if the variance matrix is diagonal, then S x = 
[l/ai.8ij] and a v = [l/a^Sij], and C xv = 0. 

A mapping T : vj <— w — (4>, J) will be linear provided 
the extent of the system in phase-space is small. Its elements 
are Tij = d-oui/dwj evaluated at (zu). Such a mapping will 
preserve the form of the distribution function. This will now 
be a Gaussian in action-angle space, with variance matrix 



&w,a 



The dynamical evolution of the system in action-angle 
coordinates is given by 



4> = 0o + n(J)f, J = constant. 

We may express 



(3) 



A w = 0(t)A WlO since A^ ,i ~ — A^j — ——AJjt, 



dJi 



and where 

e(t) = 



t 3 -n't 
o i 3 



(4) 



I3 here is the identity matrix in 3-D, and ft' represents a 
3x3 matrix whose elements are dQi /d Jj . 

Therefore the distribution function at time t is 



/(w, t) = fa exp 



where a m is now a function of time 



(5) 



(6) 



Finally, using Eq. ([5]) we may derive the distribution 
function in configuration and velocity space at time t. To 
this end, we perform a local transformation using the matrix 
T. Since this is done locally, our distribution function is still 
a multivariate Gaussian. The variance matrix at time t is 



<r ro (f) = (To©(f)T- 1 ) t a w ,o(T o 0(t)T- 1 ) 



(7) 



The variance matrix contains all the information about 
the properties of the particles on initially nearby orbits. For 
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Figure 2. Example of an orbit integrated in a Plummer potential. 

example, the evolution of the velocity ellipsoid may be de- 
rived from the velocity submatrix: a v . This submatrix de- 
scribes the velocity distribution of nearby particles at time 
t. The spatial density at a particular location x at time t 
(which is related to the spatial separation of those particles) 
is obtained by integrating the distribution function with re- 
spect to all velocities: 



p(x, t) = (2ty) foo- vl a V2 a V3 X exp 



(8) 



where a Vi=1 2 3 are the velocity dispersions along the prin- 
cipal components of the velocity ellipsoid. The matrix <r x is 
3x3, and contains all the information concerning the evolu- 
tion of the particle distribution in configuration space, in- 
cluding their separation, which is ultimately, the quantity 
that we want to measure. 



2.2 Example: Evolution in a Plummer potential 

This simple spherical gravitational potential has the form 



(r) 



GM 



Vr 2 + b 2 ' 



(9) 



Units were so chosen that G = M = b = 1 and the internal 
energy of the system is E = — |f . We define the crossing 
time of the system t cr = R/V where 7? = —GM 2 /2E and 
V 2 — —2E/M. For example, for a dwarf galaxy size system 
with b = 0.5 kpc and M = 1O 7 M , then t cr ~ 0.33 Gyr. 

We assume that the initial 6D variance matrix crj^ is 
diagonal (see Eq. [5J, with S Xj o = [l/a 2 5ij] and ov.o = 
[1/alSij], and where a x = 10 _J and a v — 10~°. We set the 
central particle of the system on an orbit whose apocentre 
is located at r a — 1.6356, as shown in Fig. [5] 

In Figure [3] we plot the evolution of the velocity dis- 
persions, the spatial density and the dispersions in config- 
uration space as function of time, for the orbit shown in 
Figure [2] These are computed using the procedure outlined 
in the previous section. 
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Figure 3. Time evolution of the velocity dispersions (top three panels on the left), spatial density (bottom left panel) and dispersions 
in configuration space (top three panels on the right), for a system moving in a Plummer potential on the orbit shown in Figure [2] The 
periodicity observed is related to the radial (and angular) orbital oscillations, as shown in the bottom panel on the right. 



Figure [3] shows that in the case of spherical potentials, 
only two of the velocity dispersions decrease in time, while 
the third one remains on average constant (it corresponds to 
the direction perpendicular to the plane of motion). These 
results imply that the configuration-space dispersions will 
increase in time, as a consequence of Liouville's theorem 
(i.e. the conservation of phase-space density). This can also 
be seen from dM ~px o xl a X2 o X3 = est. 

The form of the dispersions in velocity and in config- 
uration space has been derived explicitly in the Appendix. 
There we work in a reference frame that coincides with the 
plane of motion (this is of course possible for a spherical 
potential). In this new frame only two coordinates and two 
velocities are required to specify completely the state of sys- 
tem. In this case, the spatial density 

p oc a vl a V2 = (A„j A l , 2 )~ 1/2 , (10) 

where At, denotes the eigenvalues of the velocity submatrix 
(j v , and for which the following relation holds 

2 2 

A^A^ = {cut 4, + ust 3 + a 2 t 2 + ait + «o) • (11) 

The coefficients a* depend both on location along the orbit 
as well as on the initial extent of the system in phase-space 
(see Eg. 1X7)1. The very rapid decrease in the spatial density 
of the system observed in Fig. [3] can be understood from 
Eqs. f| 10 p and |TT} . This decrease implies a rapid separation 
of the particles (and hence of their orbits). Furthermore, the 
strong enhancements in the density seen in Fig. |3]take place 
at the orbital turning points: when p r — then X V1 x A„, — > 
and hence p — > oo. 

In the Appendix (see Eq. IA12|) . we show that the 

configuration-space dispersions a xl a X2 = \ ~T~^ — "ff - Close 



inspection of Eqs. (|lip and (|A7[) . allows us to reach the fol- 
lowing conclusions: 

• For very short timescales, the term with cvo dominates. 
In this case the separation of nearby orbits as measured by 
a x purely reflects the geometry of the orbit in phase space 
(being heavily weighted by r 2 p 2 ). 

• The terms with cti and CV4 are always positive, imply- 
ing that these will induce a rapid increase in the A„, and 
hence of the dispersions in configuration space on interme- 
diate timescales. 

• The terms with u\ and 03 can either be positive or neg- 
ative, depending on location along the orbit. This (partly) 
explains the strong oscillatory behaviour observed in Fig. [3] 

• On longer timescales, only the term cut 4 is important. 
This gives rise to the secular behaviour of density which 
decreases as 1/t 2 (as found by HW), and for the dispersions 
in configuration-space to increase as t. 

2.2.1 Relating the dispersion to the separation in 
configuration space 

Our main aim is to study the separation of initially nearby 
orbits in configuration space. The question is how this sepa- 
ration is related to the dispersions or the variances (i.e. the 
inverse of the eigenvalues) in the configuration-space matrix 

We examine three possibilities obtained by performing 
three different kinds of averages of the configuration-space 
dispersions: 

(i) the geometric mean: A 9 = (a xl a X2 a X3 ) 1 ^ 3 , 

(ii) the arithmetic mean: A a = (a xi + a X2 + a X3 )/3, 

(iii) the modulus: A m = yj a xi + a\ 2 + cr X3 /3. 

Figure[4]shows the behaviour of these different averages. 
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Figure 4. Time evolution of three possible averages of the disper- 
sions in configuration-space, as obtained through our formalism. 
Note the very fast initial growth and the linear behaviour at late 
times. The dashed curve in each panel represents the average sep- 
aration (A r ) of 1000 nearby orbits. Note the excellent agreement 
between (A,-) and A a (middle panel). 

In all cases, one observes initially a very fast increase in the 
measured values, while at late times the growth proceeds 
linearly with time [see also Eqs. (|A6[) and ()A10|| ]. 

The question now is how to relate the above defined 
averages to the separations between two initially nearby or- 
bits. To address this we have generated 1000 orbits with 
initial conditions distributed according to a Gaussian in con- 
figuration space with dispersion a x — 10 -5 and a v = 10 -5 
around the orbit shown in Fig. [3] We measure the separation 
A r = ri — ro between this orbit and the 1000 neighbour- 
ing trajectories, and derive the average (A r ). This is plotted 
as a dashed curve in the panels of Fig. [3] As can be seen, 
the arithmetic mean A a of the configuration-space disper- 
sions computed using the analytic formalism discussed in 
the previous section provides an excellent measurement of 
the separation between nearby orbits. 

Figure [4] shows that the separation of nearby orbits in 
smooth integrable potentials exhibits an initial rapid diver- 
gence, which is followed by a secular increase which is linear 
in time. Note that this occurs for a completely integrable sys- 
tem, without any degree of chaos. Therefore we see already 
that the initial exponential divergence cannot be attributed 
to any form of chaos whatsoever. It simply reflects the way 
an orbit evolves in phase-space, as shown in the Appendix. 

2.2.2 Direct comparison to N-body simulations 

The behaviour visible in Figure [4] is strinkingly similar to 
th at observed in th e N-bo dy simulations shown in Fig. 1 
of IValluri fc Merrittl ifeOQOl) and in t he frozen N-body real- 
izations of lKandrup fc Siderid (120031 ) . To provide the reader 
with a more direct comparison to our analytic results, we 
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Figure 5. Time evolution of the average separation (A r ) of 100 
nearby orbits integrated in a frozen N-body realization of the 
Plummer sphere (solid curve). The error bars correspond to the 
error on this average. The dashed curve represents the separation 
of nearby orbits as measured by A a using our analytic prescrip- 
tion, while the dotted line represents the same quantity but esti- 
mated from the 100 orbits integrated in the N-body realization. 

will now analyse the divergence of an ensemble of initially 
nearby orbits evolved in a potential represented by a finite 
number of particles. 

We have generated a realization of the Plummer sphere, 
whose density profile is 

which we have truncated at a radius rt = 12.1976, that 
encloses 99% of its mass. We represent this system with 
N — 128, 000 par ticles, and use a numeric al softening 
e = 0.025& (see e.g . lAthanassoula et all (|200Ch ). Following 
iKandrup fc Siderid l|200ll ) we integrate orbits in this frozen 
(in time and space) N-body realization. The integration of 
the orbits was performed using a Runge-Kutta-Fehlberg al- 
gorithm of order 4-5. 

As in the previous section, we follow the evolution of 100 
orbits distributed according to a Gaussian in phase space 
with initial dispersion a x — 10~ 5 and a v — 10~ 5 around the 
orbit shown in Fig. 2. For this ensemble we have measured 
the time evolution of the average separation (A r ) and of the 
6D variance matrix. 

The results are shown in Fig. [S] The initial behaviour is 
very similar, whether derived using our analytic formalism 
(A a , dashed curve) or using the average separation of 100 
orbits in the N-body representation of the system ((A,-), 
solid curve). The small differences can be attributed to two 
causes. First of all, the finite sampling of phase-space around 
the central orbit introduces an error in the average, which is 
quantified by the error bars shown in this figure. Secondly, 
(Ar) is not exactly identical to A a (see Figure |4}. If we 
compare the quantity A a obtained using the 6D variance 
matrix from the frozen N-body simulation (dotted curve) 
with that from the analytic formalism (dashed curve), this 
difference almost disappears. The two curves are virtually 
indistiguishable over a timescale of a few crossing times. 

2.2.3 Dependence on initial conditions 

It is also interesting to understand how the separation of 
initially nearby orbits depends on their initial conditions. In 
particular, how the initial divergence depends on the region 
of phase-space sampled at short times. 
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Figure 6. Time evolution of the separation of nearby orbits as 
measured by the arithmetic mean A a for the same orbit discussed 
in previous figures. The sensitivity to initial location along a given 
orbit is evidenced by the various curves: solid corresponds to ini- 
tial location at pericentre; short-dashed to apocentre and long- 
dashed to the average distance between these turning points. 

In Figure [6] we plot the time evolution of the arithmetic 
mean of the three dispersions in configuration space A a ob- 
tained using our analytic prescription for the orbit shown 
in Fig. [2] We now plot the behaviour for different starting 
points along this orbit: apocentre, pericentre and (apocentre 
+ pericentre)/2. We see clearly that the behaviour at short 
times depends on the initial location along the orbit. The 
initial divergence is in all cases nearly exponential, but has 
largest amplitude (it lasts longer) when the integration is 
started near pericentre. 

It is also interesting to study the behaviour of different 
sets of nearby orbits. Figure [7] shows the evolution for two 
additional examples. The top panel corresponds to an or- 
bit constrained to move in the inner regions of the system 
(pericentre r v = 0.626 and apocentre r a = 0.886), while the 
bottom panel has r p — 1.736 and r a = 2.946, and hence it 
is constrained to the outskirts. Clearly the amplitude of the 
initial growth phase depends on the regions of phase-space 
the orbits probe. 

Note that, because the quantities shown are normal- 
ized to their initial conditions, these results are independent 
of the initial separation of the orbits (or the values of an 
in our formalism). This is perhaps, the most characteristic 
difference between an integrable and a chaotic system. The 
amplitude (or the rate) of the initial divergence for a given 
orbit is always the same in the integrable case, irrespective 
of initial separation. To the contrary, in a chaotic system, 
in the limit of infinitesimal perturbations, the orbits may be 
trapped near a resonance and cease to be chaotic to become 
regular. Hence the amplitude of the initial divergence will, 
in the chaotic case, depend strongly on the initial separation 
of the orbits. 



2.2-4 Miller's instability and the initial behaviour 

The above analysis shows that the initial very rapid diver- 
gence of nearby orbits is a generic feature of dynamical sys- 
tems. It is not only observed in the Plummer potential dis- 
cussed here, but also in all integrable potentials studied by 
HW (e.g. Fig. 7 and 9 of their paper). 

The initial behaviour is nearly exponential, as shown 
in Fig. [7] for the orbits discussed so far. The rate of diver- 
gence -the equivalent of the "short-term" Lyapunov expo- 
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Figure 7. Time evolution of the separation of nearby orbits as 
measured by the arithmetic mean A . The various panels repre- 
sent orbits probing different regions of the system, all integrated 
from their apocentres. The dotted curves are exponential fits to 
the initial (transient) behaviour, while the dashed curves are lin- 
ear fits to the long-term (secular) behaviour of A a . 

nent, is Xe ~ 13.6/t cr for the inner orbit, Xe ~ lQ.7/t cr and 
Xe ~ 17.57/t C r for the intermediate and outer orbits, respec- 
tively. In all cases, the secular behaviour at late times can 
be fit by a linear function of time A a = A aj o + t/t sec , where 
the divergence timescale is t sec ~ 1.3t cr , t sec ~ 0.77t cr and 
t sec ~ Alter for the different orbits. 



3 DISCUSSION 

Our analysis shows that the initial nearly-exponential di- 
vergence of nearby orbits in N-body systems is not due to 
chaos. It is present also in integrable smooth potentials, and 
it reflects a power-law divergence modulated by the shape 
of an orbit in phase-space. 

It is interesting to note that the rates of divergence that 
we measure using our form alism are in very good agree- 
ment with those obtained by Irlemsendorf fc Merrifra (|2002T ) 
for N-body realizations of the same Plummer sphere. These 
authors find a characteristic e-folding time of t cr /20 for sys- 
tems with N ~ 10 J particles, which is very comparable to 
the values obtained in Section [2.2.41 They also find a weak 
dependence on N, which may also be readily understood 
within our framework. Such a dependence is induced by the 
very rapid decrease in the spatial density of the system. If 
a "relatively" small number of particles is used in a N-body 
simulation, then the density cannot be mapped properly. For 
example, to measure a decline in the density of on a 

timescale of ~ 3t cr as observed in Fig. [3] N-body realizations 
with at least 10 nearby particles are neede d. 

Previous wor ks, including I Miller! (| 19641 ) and 

iKandrup fc Siderisl (|2003l ) have also noted an oscilla- 
tory behaviour in the divergence of nearby orbits. Our 
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analysis, as well as Figure [3] show that this is due to the 
modulation produced by the periodicity of a r e gular orbit 
in phase-space. It is not, as suggested bv lMillerl (jl964n . due 
to the formation of tight binaries in an N-body system. The 
fact that such behaviour was visible in the various N-body 
studies presented in the literature, in fact demonstrates 
that such N-body systems were faithful representations of 
the true (int egrable) system , at least on short timescales . 
As stated bv lKandrupl dHH) and lValluri fc MerrittJ (|2000h . 
the Lyapunov exponents need to be measured in the 
limit of infinite time intervals; short-time exponential-like 
divergences do not imply chaotic behaviour. 
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APPENDIX A: MATRICES IN A SPHERICAL POTENTIAL 

For spherical potentials <&(r), we may choose a system of coordinates that coincides with the plane of motion of the system. 
In this plane the position of a particle is specified by its angular (ip) and radial (r) coordinates. The actions of an orbit in 
this case are: 

r r 2 

(Al) 



L = J^=p^, J r = - dr- y / 2[E - $( 



■)]» 



L 2 , 



where L is the total angular momentum of the particle, E is its energy, and ri and the orbital turning points. 

In order to track the evolution of the dispersions of our initial distribution function, f(vj,to), we perform the following 
sequence of operations. Firstly, we transform from Cartesian coordinates zu = (x, v) to action angle variables w = (6,3). The 
distribution function is then evolved in this space, after which, we transform back to Cartesian coordinates (see Figure [IJ. 

For the sake of simplicity, here we begin with a distribution function already expressed in terms of the action-angle 
variables and we also assume that, initially, the variance matrix is diagonal, i.e., a W! o = [ou5ij\. With the time evolution 
operator, &(t), known, we can compute the variance matrix at any given time t as a w (t) = ©(*)' cr W! o@(t). Equation 
shows @(t) for the 3-D case, but we reduce the equation for our purposes to the 2-D case. 

After evolving the system in the action-angle space we need to transform back locally to configuration and momenta 
space lj — (x, p) using the transformation matrix T _1 . The elements of this matrix are related to the second derivatives of 
the characteristic function W(q,3). In our case 



(A2) 



1 


*12 


tl3 


tl4 





*22 


*23 


*24 








1 








*42 


*43 


*44 



with 



h(r), Tr K 

tl2= ±- L W 3 A + —, *13 

h(r) Q r 

t22 = — VV4A H , *23 

h(r) 

*42 = ^ — , *43 



where 



W33 + ^34*43, 
K 



tl4 = W / 34*44, 

t24 = W^nt 44, 
, Pr 



h(r) = -&(r) + —, Pr 
and 



2[E-Mr)]-^, 



fit* r, 



W33 



Wa. 



w 34 = 



d 2 W 
OL 2 

d 2 W 

dJr 2 = 

d 2 W 
OLdJr 



r dr ( 9fL 



i p r \ dJ^p r 2 
dr ( d0, r 

Pr \ dJ r Pr 



dr f dQ,^ k 

Pr \ 8Jr Pr 



n 2 ' 

Pr 



Subindices 1 and 3 in the expressions above refer to directions associated with if), such as, 
are related to r. For more details about this procedure we refer the reader to HW. 
Given T _1 , the variance matrix at time t is expressed as: 

where the elements t%j are evaluated at (x(t)). Substituting T _1 , ®(t) and <rjj, in the above expression, then 



and J^, whereas 2 and 4 



o a (t) = 



(711 

{1,2} 
{1,3} 
{1,4} 



criiyl 

(J22D 2 + (744*42 

{2,3} 
{2,4} 



011B 

anAB + (T22DE + Gat 42*43 



G\\AC + G22DF + a"44t42*44 



anB + 022E + (T33 + (744*43 a\\BC + 022EF + 044*43*44 



{3,4} 



&11C + U22F + a 44* 44 



(A3) 



(A4) 



where 
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A = ti2 — ^34t42i, B = tn — (fl 33 — r2 34 i43)f, 
D =t22 ~ Q! 44 t 4 2~t, E = t23 — (^34 — ^ 4 4t43)f, 



C — tl4 — fi^f^t, 
-F = t24 — f^44t44i- 



In general, one is more interested in the properties of the debris in velocity space, rather than in momenta space. Therefore 
we transform the variance matrix according to cr ro (t) = Tp^ v o-a(t)T p —, v , with 



1 

1 

Wy, 









r 

1 



To obtain an expression for the time evolution of the velocity dispersions we focus our attention on what happens around 
a particular point (x(t)) in configuration space located on the mean orbit of the system. This is equivalent to studying the 
velocity submatrix of the variance matrix <r OT (t), that is 



r 2 (auB 2 + 022E 2 + 033 + 0-44*43) r{a xl BC + a 22 EF + 044*43) 



{1,2} 



OliC + (J22F + {744*44 



(A5) 



By diagonalizing the matrix ov we obtain the principal axes of the velocity ellipsoid at the point (x(t)), and the associated 
dispersions. The eigenvalues of o v are the roots of the characteristic equation: det[o v — AT] = 0. An interesting quantity is 
for example, A V1 A V2 because it is inversely proportional to the density: p tx a vl a V2 — (X V1 A„ 2 ) _ • In our case: 



A^i A^2 
where 

Q4 = 

Q 3 = 
Q.2 = 

Ql = 
OLQ = 



2 2 

r p r 



(a?4* 4 + a 3 t 3 + a2t 2 + ait + cto) , 



(A6) 



o"n(J22(det fi') 2 , 

2aiiCT22 det f2' (2W 3 4« 3 4 - ^33^44 - W i4 U i3 ), 

ana 2 2 (2 det fl' det W + {Q' i4 W 33 + Q' 33 Wm) 2 + 4Wm(^mWs* ~ ft^Q'^Wu - 0^044 W33)) + 

(o"H(J33 + O22O44) ^'34 + <Tll<744f2' 33 + 0"22CT33fi'44, 

2(0-110-22 det W (2f2 34 VK34 - fi 44 VK33 - O33W44) - ^34^34 (0110-33 + 0220-44) - 01 1044 033^33 - 
022033SI44W44) , 

(on022)(det W) 2 + W34 (011O33 + O-22O44) + 0-nO44H / 33 + O22O33W44 + 0"33O44, 



(A7) 



with 

detW = W33W44 - W% A . 

These equations explicitly show the behaviour of principal axes velocity dispersions: 

• For very short timescales, the term with ao dominates. In this case the behaviour purely reflects the geometry of the 
orbit in phase space (being heavily weighted by r 2 p 2 ). 

• The terms with 02 and 04 are always positive, implying that these will induce a rapid increase in the A„, or a rapid 
decrease of the velocity dispersions on intermediate timescales. 

• The terms with ai and 03 can either be positive or negative, depending on location along the orbit (i.e. the Wij vary in 
magnitude and sign). This explains the strong oscillatory behaviour observed in Fig. [3] 

• On longer timescales, only the term a 4 t 4 is important. This gives rise to the secular behaviour of density which decreases 
as 1/t 2 , and the velocity dispersions to behave as 1/t for long timescales. 

To obtain the expression for the time evolution of the dispersions in configuration space we integrate the distribution 
function with respect to all velocities (see Eq. |SJ). In practice, we first transform cr w (t) from polar to Cartesian coordinates, 
oi(i) = (TOW^T', where 



T' 



sin(f/j) 
r 

cos(V') 
sin(^)p r 

r 



cos(?/>) 
r 

COs(lp)p r 

r 

cos(ip)v^ 






— sin(tp) cos( , 0) 

cos(V0 sin(f/>) 



(A8) 



We express o^ as 
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A 

B + 



B 

C 



where the 2x2 matrices A, C and B represent the position submatrix, the velocity submatrix, and the cross correlation 
between positions and velocities, respectively (as in Eq. |SJ. Then, the matrix cr x is obtained from the integration of the 
distribution function over the velocities: 



511 S12 

512 S22 

where the elements Sij are related to the dispersions in configuration space. These elements can be expressed as: 

_ det Tjj 
StJ ~ detC ' 
with 



(A9) 



dij 


bu 


b,2 




Cll 


C12 


bj2 


Cl2 


C22 



where a%j, &y and cy are elements of the matrices A, B and C respectively. The diagonalization of the matrix <j x yields the 
values of the dispersions along the principal axes of the system in configuration space since <j Xi = 1/a/ A r i , where \ Ti are the 
eigenvalues of a x . 

Solving the characteristic equation for cr x we finally obtain: 



where 

ft = 



(2A„ 1 A l , 2 )~ 1 j3 2 t 2 + Pit + ft ± y/ifot 2 + Pit + /3 ) 2 - 4X V1 A„ 2 det ag 



(A10) 



0\\022T 2 [((T44^'34 + 4\) (p 2 + r 2 — 2^34(044^33 + (T3 3 Q.' i4 )r 2 KQ r + (044^33 + 34 )r 2 Q.l] , 

01 = -2on(J22r 2 [((J44H / 34fi34 + ^33 W44Q44) (Pr + r 2 *, 2 ) + (((733W44 + 044^33)^34 + W^O^^M + ^33^44 ))r 2 |cfi r 
-((744 1^33^33 + ^33 W 34^34) r 2 ^lV\ , 

00 = r 2 [an (<T220-3 3 W44 + 044 (022 W34 + 033)) (p 2 r + r 2 k 2 ) - 2<7n (T22W34 (0-33 Wu + 044 W^r 2 'K,Sl r + 

(T22 (o-ii(T33 W34 + <T4a{cti\W% 3 + 033 )) r 2 ^ 2 .] . (All) 
Finally, multiplying both eigenvalues: 
det a?„ 



A-u^ A^ 2 



(A12) 



